--- title: "Cholesky Factor Investigation" permalink: /Research/CholeskyFactorInvestigation/ author_profile: true ---
In this notebook I look at a few (4) cholesky factors, used to whiten KSE data prior to implementing a least dquares solver to compute the Wiener filter.
The sizes of the Cholesky factors are 50, 100, 200, and 500.
Below is a section for each size. Within each section are a number of plots of (1) the the powerspectrum of various columns of the whitened design matrix and (2) the above diagonal columns of the Cholesky factors (C, C100, C200, C500). Some section also contain calculations for the infinity norm of the difference between corresponding submatrices of the various Cholesky factors.
The data suggests that
using JLD
using LinearAlgebra
using Statistics
using PyPlot
at = include("../../../../Tools/AnalysisToolbox.jl")
function get_hash(pred::Array{<:Number,2};M_out = 30) where T<: Number;
nu, steps = size(pred)
pred = Array(transpose(pred))
PRED = zeros(ComplexF64,steps-M_out+1,M_out*nu)
for m = 1:M_out
PRED[:,nu*(M_out-m)+1:nu*(M_out-m+1)] = @view pred[m:(steps + m - M_out),:]
end
PRED = reverse(PRED,dims = 2)
PRED
end
# Load Data
#This is writen to allow for loading any of 5 different runs of the KSE solution
geni = 1
gen = "lin1e5_r$geni" # this is just a reference designation it shows up in the
# output file. I think of generatrion.
savefile_pred = "../../../../../data/KSE_Data/ks_pred_$gen.jld"
println("Sol load location: " * savefile_pred)
Pred = load(savefile_pred,"pred");
Sig = load(savefile_pred,"sig")
sig = Sig[1:1,:]
pred = Pred[1:1,:]
size(pred), size(sig)
function get_cholfact(pred, M_out; pcond = 1, eps = 1e-5, verb = true)
Xh = get_hash(pred; M_out = M)
w = I
for i = 1:pcond
σ = eps^2*var(Xh[:,M])
A = Xh'Xh/size(Xh,1) + σ*I
verb && println("preconditioning $i original condition number: ",cond(A))
println("Information on the factoring of the approximated covariance matrix:")
@timev L = cholesky(A).L
verb && println("preconditioning $i Cholesky condition number: ",cond(L))
W = inv(conj(L))
w = W*w
Xh = Xh*transpose(W)
verb && println("preconditioning $i new condition number: ",cond(Xh))
end
C = transpose(w);
Xh, C
end
M = 50
Xh, C = get_cholfact(pred, M);
function plot_proc(Xh, M_out; Nex = 2^13, Θ = 2π*(0:Nex-1)/Nex)
for j = 1:10:M_out
figure()
for i = 0:3:10
title("Power spectrum of various columns of whitened design matrix ($M_out × $M_out)")
S = at.z_crossspect_dm(Xh[:,i+j],Xh[:,i+j]; Nex)
loglog(Θ,real(S), label = "Column $(i+j)")
ylabel("Power")
xlabel("Frequency")
legend()
end
end
end
plot_proc(Xh, M)
function plot_WF(C, M_out)
for j = 1:10:M_out
figure()
for i = 0:3:10
title("Relevent portions (diagonal upwards) of various columns Cholesky factor ($M_out × $M_out)")
plot(real(C[i+j:-1:1,i+j]), label = "Column $(i+j)")
ylabel("Real part of coefficient")
xlabel("lag")
legend()
end
end
end
plot_WF(C, M)
M = 100
Xh, C100 = get_cholfact(pred, M);
plot_proc(Xh, M)
plot_WF(C100, M)
function plot_WF_comp(C1,C2)
m1 = size(C1,1)
m2 = size(C2,1)
for j = 1:10:min(m1,m2)
figure()
for i = 0:3:10
title("Overlay of whitening filters from corresponding columns of ($m1 ×$m1 solid) and ($m2 × $m2 dotted) cholesky factors")
plot(real(C1[i+j:-1:1,i+j]), label = "Column $(i+j)")
plot(real(C2[i+j:-1:1,i+j]),":", label = "Column $(i+j)")
ylabel("Real part of coefficient")
xlabel("lag")
legend()
end
end
end
plot_WF_comp(C,C100)
M = 200
Xh, C200 = get_cholfact(pred, M);
plot_proc(Xh, M)
plot_WF(C200, M)
plot_WF_comp(C,C200)
plot_WF_comp(C100,C200)
norm(C - C100[1:50,1:50], Inf)
norm(C - C200[1:50,1:50], Inf)
norm(C100 - C200[1:100,1:100], Inf)
norm(C100[1:50,1:50] - C200[1:50,1:50], Inf)
norm(C100[51:100,51:100] - C200[51:100,51:100], Inf)
norm(C200[101:150,101:150] - C200[151:200,151:200], Inf)
var(C[:]),var(C100[:]), var(C200[:])
M = 500
Xh, C500 = get_cholfact(pred, M);
plot_proc(Xh, M)
plot_WF(C500, M)
plot_WF_comp(C,C500)
plot_WF_comp(C100,C500)
plot_WF_comp(C200,C500)
var(C[:]),var(C100[:]), var(C200[:]), var(C500[:])
var(real(C[:])),var(real(C100[:])), var(real(C200[:])), var(real(C500[:]))
var(imag(C[:])),var(imag(C100[:])), var(imag(C200[:])), var(imag(C500[:]))
mean(imag(C[:])),mean(imag(C100[:])), mean(imag(C200[:])), mean(imag(C500[:]))